setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/20210103.motif.analysis")
file=read.csv("807.psy.ASH.add.eQTL.GWAS.add.DEG.LIBDeQTL.HMM.state.motif.csv",header=T)
file1=file[file$BF_in_DC>10,]
file1=file1[file1$eqtl.no.na.num>0|!is.na(file1$LIBD.eQTL.type),]
file1=file1[file1$gwas.no.na.num>0,]

tf=read.csv("cor.mean.vaf.vs.score.allele.diff.csv",header = T)
id=tf[tf$cor.p.value<0.05&abs(tf$cor.estimate)>0.7,]$TF

tmp=file1[file1$geneSymbol %in% id,]

#以下是提取4个位点的VAF
setwd("E:\\5hmc_file\\2_5hmc_yjp_bam\\ASM\\bayes_pvalue_beta0")
group1=c("X2B_X1T","M8_M7","M6_M5","M2_M1","M48_M47","M50_M49","M28_M27","M30_M29","M26_M25","M35_M36","M18_M17","M20_M19","M22_M21","M40_M39")

id=tmp$unitID
i=1
fn1=paste0(group1[i],".bayes_p.txt")
file=read.table(fn1,head=T,sep = "\t")
file$unitID=paste(file$chrom,file$position,file$ref,file$var,sep=":")
file=file[file$unitID %in% id,]
file=data.frame(unitID=file$unitID,file$normal_reads1,file$normal_reads2,file$normal_var_freq,normal_bayes_beta0=file$normal_bayes_beta0,normal_bayes_pvalue=file$normal_bayes_pvalue,file$tumor_reads1,file$tumor_reads2,file$tumor_var_freq,tumor_bayes_beta0=file$tumor_bayes_beta0,tumor_bayes_pvalue=file$tumor_bayes_pvalue)
file$file.normal_var_freq=round(as.numeric(gsub("%","",file$file.normal_var_freq))/100,digits = 4)
file$file.tumor_var_freq=round(as.numeric(gsub("%","",file$file.tumor_var_freq))/100,digits = 4)
str1=unlist(strsplit(group1[i],"_"))
str2=c(".reads1",".reads2",".var_freq",".bayes_beta0",".bayes_pvalue")
names(file)=c("unitID",paste(str1[1],str2,sep=""),paste(str1[2],str2,sep=""))
rt=file
for(i in 2:14){
  fn1=paste0(group1[i],".bayes_p.txt")
  file=read.table(fn1,head=T,sep = "\t")
  file$unitID=paste(file$chrom,file$position,file$ref,file$var,sep=":")
  file=file[file$unitID %in% id,]
  file=data.frame(unitID=file$unitID,file$normal_reads1,file$normal_reads2,file$normal_var_freq,normal_bayes_beta0=file$normal_bayes_beta0,normal_bayes_pvalue=file$normal_bayes_pvalue,file$tumor_reads1,file$tumor_reads2,file$tumor_var_freq,tumor_bayes_beta0=file$tumor_bayes_beta0,tumor_bayes_pvalue=file$tumor_bayes_pvalue)
  file$file.normal_var_freq=round(as.numeric(gsub("%","",file$file.normal_var_freq))/100,digits = 4)
  file$file.tumor_var_freq=round(as.numeric(gsub("%","",file$file.tumor_var_freq))/100,digits = 4)
  str1=unlist(strsplit(group1[i],"_"))
  str2=c(".reads1",".reads2",".var_freq",".bayes_beta0",".bayes_pvalue")
  names(file)=c("unitID",paste(str1[1],str2,sep=""),paste(str1[2],str2,sep=""))
  rt=merge(rt,file,by="unitID",all=T)
}
write.csv(rt,"../20201216/4.Coordinates.site.csv",quote=F,row.names = F)
